//@+leo-ver=5-thin
//@+node:sonam.20180318222433.1: * @file FivePlanet.pas
//@+at
//  該單元實現了五曜的計算
// 
//@@c
unit FivePlanet;

//@+others
//@+node:sonam.20180328224215.1: ** Interface
interface

uses SysUtils, FracData, PlanetCycle, SysElmnt;

Type
//@+node:sonam.20180328224249.1: *3* 類聲明
	EZeroYearException = class(Exception);
	TFivePlanet = class(TObject)
    public
		constructor Create(bk: CalcBook);
        function CommonDay() : Integer;               //公日
        function RespectiveDay(p: Planet) : Integer;  //殊日
        function MeanRevolution(p: Planet) : TPlanetData; //武遲文步
        function MeanParallax() : TSunPlanetData;         //文遲武步
        function SlowStep(p : Planet) : TCompFraction; //遲步
        function SlowMove(p : Planet) : TCompFraction; //遲行定數
        function FastStep(p : Planet; var SecondHalf,
                            backStep : Boolean) : TGenPlanetData;   //疾步
        function FastMove(p: Planet; var Direction : Char) : TGenPlanetData;  //疾行定數

    private
		fBook : CalcBook;

		fYear : Integer;
		fMonth : Integer;
		fDay : Integer;
        fLeapMonth : Boolean;
		fCycle : TPlanetCycle;
		procedure SetYear(y: Integer);
		procedure SetMonth(m : Integer);
		procedure SetDay(d: Integer);
        procedure SetLeapMonth(l: Boolean);
        function AccumSunDays(p : Planet) : Int64;
        function GetMeanSlowMove(p : Planet) : TCompFraction;
        function PlanetAccumMonth(y : Integer) : Int64;
	public
        property Year : Integer read fYear write SetYear;
        property Month: Integer read fMonth write SetMonth;
		property Day : Integer read fDay write SetDay;
        property LeapMonth : Boolean read fLeapMonth write SetLeapMonth;

        //此處的y雖然有勝生年、西元紀年和自丁卯已過週年，m和d都是藏歷月和日
		procedure SetDate(y, m, d: Integer);
        procedure SetChristDate(y, m, d: Integer);
        procedure SetPastDate(y, m, d: Integer);

		property Book : CalcBook read fBook;
		property cycle : TPlanetCycle read fCycle;
	end;
//@+node:sonam.20180328224405.1: ** Implementation
Implementation

Const
//@+node:sonam.20180318221359.1: *3* 各曜的近地點
//@+at 各曜的遠地點
// 
//@@c
    Apogee : array [0..4, 0..5] of Integer =
        (   (9,30,0,0,0,0),
            (12,0,0,0,0,0),
            (18,0,0,0,0,0),
            (16,30,0,0,0,0),
            (6,0,0,0,0,0)  );
//@+node:sonam.20180328224744.1: *3* 基本方法
constructor TFivePlanet.Create(bk: CalcBook);
begin
	fBook := bk;
    fCycle := TPlanetCycle.Create();  
    LeapMonth := false;  
end;

procedure TFivePlanet.SetYear(y: Integer);
begin
	fYear := y;
end;

procedure TFivePlanet.SetMonth(m: Integer);
begin
	fMonth := m;
end;

procedure TFivePlanet.SetDay(d: Integer);
begin
	fDay := d;
end;

procedure TFivePlanet.SetLeapMonth(l: Boolean);
begin
	fLeapMonth := l;
end;

procedure TFivePlanet.SetDate(y, m, d: Integer);
begin
	fYear := y;
	fMonth := m;
	fDay := d;
end;

procedure TFivePlanet.SetChristDate(y, m, d: Integer);
begin
    if y = 0 then
      raise EZeroYearException.Create('Zero year');

    if y > 0 then
      fYear := y + ChristianBeginYear
	else
      fYear := y + ChristianBeginYear + 1;

	fMonth := m;
	fDay := d;
end;

procedure TFivePlanet.SetPastDate(y, m, d: Integer);
begin
	fYear := y + ManiBeginYear;
	fMonth := m;
	fDay := d;
end;

function TFivePlanet.PlanetAccumMonth(y : Integer) : Int64;
begin
    result := PlanetCycle.PlanetAccumMonth(y, month, LeapMonth);
end;
//@+node:sonam.20180318221217.1: *3* 公日
//@+at 公日即是所求日距歷元的太陽積日，五曜公用，故名公積日。
// 這裏因為採用三日共週年計算，后面的步驟中不會用到。
// 
//     <公日> :=
// 
//@@c    
function TFivePlanet.CommonDay() : Integer; // 殊日
begin
    result := 0;
end;
//@+node:sonam.20180318221233.1: *3* 殊日
//@+at 
// 殊日即五曜經過白羊宮首后已經運行的日數，以五曜各不相同，故名“殊日”。殊是相对
// 五曜所共同的公日而言的，公日即是輪值曜的太陽積日。此處採用《總滙如意寶》140頁到141
// 頁所說的計算方法。
// 
//     <各曜太陽積日> := (<各曜太隂積日>) * 11135 / 11312
// 
// 這裏的11135/11312是太隂日化為太陽日的因子，此處要在有餘數時加一。而其中:
// 
//     <各曜太隂積日> := (<各曜三日共週年> * 12 + <已過月>) * 67 / 65 * 30 + <日期>
// 
// 由太陽積日就可以得到殊日：
// 
//     <各曜殊日> := <各曜太陽積日> mod <各曜的太陽日週期>
// 
//@@c
function TFivePlanet.AccumSunDays(p: Planet) : Int64;
Var
    pastLunarDays : Int64;
    pastDays : Int64;
begin
	pastLunarDays := PlanetAccumMonth(cycle.YearInCycle(book, p, year)) * 30 + day;
    if pastLunarDays * 11135 mod 11312 = 0 then
        pastDays := pastLunarDays * 11135 div 11312
    else
        pastDays := pastLunarDays * 11135 div 11312 + 1;
    result := pastDays;
end;

function TFivePlanet.RespectiveDay(p: Planet) : Integer; // 殊日
begin
    result := AccumSunDays(p) mod cycle.DayCycle(p);
end;
//@+node:sonam.20180318221245.1: *3* 武遲文步
//@+at
// 武遲文步，是武曜的遲行中數和文曜的引檢步數的簡稱。它在天文學上的意義就是各曜的
// 平均本身行。
// 
// 對於武曜，由於其公轉週期比地球的公轉週期慢，地球運行引起的視差移動相對於曜的本
// 身行所引起的視差移動要快，因此將武曜的本身行所引起的視差移動稱為遲行。此處計算
// 的是尚未考慮遲步的中間結果，因此稱為“遲行中數”
// 
// 武曜的遲行中數和文曜的引檢步數，即所求日該曜距春分點的弧長。引檢步數減去遲行中
// 數后，爲查檢遲行步度盈縮表時用的引數，故名"檢步"。
// 
// 根據《總滙如意寶》第142頁，武遲文步的計算方法是：
// 
//     <武遲文步> := <公積日> * <各曜的平均運行速度>
// 
//@@c
function TFivePlanet.MeanRevolution(p: Planet) : TPlanetData; // 武遲文步
Const
//@+at
// 各曜的平均太陽日運行速度表
//@@c
    PlanetVelocityTab : array [0..4, 0..4] of Integer =
        (   (0,2,21,2,208),
            (0,0,22,2,226),
            (0,0,9,0,918),
            (0,18,24,5,4687),
            (0,7,12,3,345)  );
begin
    result := TPlanetData.Create(p, PlanetVelocityTab[Ord(p)-Ord(Mars)]);
    result.Multiply(RespectiveDay(p));
    result.Calibrate();
end;
//@+node:sonam.20180318221257.1: *3* 文遲武步
//@+at 所謂的文遲武步，即文曜的遲行中數和武曜的檢步序數的簡稱。它在天文學上的
// 意義即地球的公轉引起的平均視差移動，也即計算五要素的過程中算出的中日。
// 
// 這里採用《總滙如意寶》第142頁到143的計算方法：
// 
//     <文遲武步> := <日之太陽積日> * TSunPlanetData(0,4,26,0,41,1849)
// 
// 而其中：
// 
//     <日之太陽積日> := <公日> mod <日之太陽日週期>
// 
//@@c
function TFivePlanet.MeanParallax() : TSunPlanetData;
begin
    result := TSunPlanetData.Create([0,4,26,0,41,1849]);
    result.Multiply(AccumSunDays(Week) mod cycle.DayCycle(Sun));
    result.Calibrate();
end;
//@+node:sonam.20180318221416.1: *3* 遲行步度表
//@@c
function GetSlowStepFactor(p : Planet; zodiac : Integer) : Integer;
Const
    SlowStepFactorTab : array [0..4, 0..5] of Integer =
        (   (-25, 25, 18, 7, -7, -18),
            (-11, 11, 9, 3, -3, -9),
            (-22, 22, 15, 6, -6, -15),
            (-10, 10, 7, 3, -3, -7),
            (-5, 5, 4, 1, -1, -4)   );
Var
    i : Integer;
begin
    if zodiac = 5 then
        i := 0
    else
        i := zodiac + 1; 
    result := SlowStepFactorTab[Ord(p)-Ord(Mars)][i];
end;

function GetSlowStepAccum(p : Planet; zodiac : Integer) : Integer;
Const
    SlowStepAccumTab : array [0..4, 0..5] of Integer =
        (   (0, 25, 43, 50, 43, 25),
            (0, 11, 20, 23, 20, 11),
            (0, 22, 37, 43, 37, 22),
            (0, 10, 17, 20, 17, 10),
            (0, 5, 9, 10, 9, 5) );
begin
    result := SlowStepAccumTab[Ord(p)-Ord(Mars)][zodiac];
end;
//@+node:sonam.20180318221432.1: *3* 遲步
//@+at 
//  所謂遲步，是各曜的遲行定數与中數之間的快慢脩正值，其計算方法在《總滙如意寶》
// 128至129頁：
// 
//     <遲步> := <遲行積步> + <各曜凈行>
// 
// 用遲行中數減去遠地點黃經之后，化為弧刻，然后除以每宮運行的弧刻數135，得到的
// 商是宮序，餘數是該曜距所在宮宮首的弧刻數。以宮序作為下標查檢遲步表，得到盈縮
// 積和損益率。
// 
//     <各曜凈行> := <損益率> * TPlanetData(0,餘數,遲行中數的剩餘位)
// 
// 由於武曜使用的是五位的分數，而文曜使用的是六位的分數，這裏將使用其基類TCompFraction
// 作為函數返囬的數據類型。
// 
//@@c
// 將曜的遲行中數減去其遠地點，得到從遠地點開始一個週期內已運行的弧長
function ArcLenFromApogee(p: Planet; pd: TCompFraction) : Boolean;
Var
    HalfCycle : TCompFraction;
begin
    pd.Sub( CreateOuterMoveData(p, Apogee[Ord(p)-Ord(Mars)]) );
    HalfCycle := CreateOuterMoveData(p, [13,30,0,0,0,0]);
    if pd.GreatEqual( HalfCycle ) then
    begin
        pd.Sub( HalfCycle );
        result := true;
    end
    else
        result := false;
end;

function OuterPlanet(p : Planet) : Boolean;
begin
    result := (p = Mars) or (p = Jupiter) or (p = Saturn);
end;

function TFivePlanet.GetMeanSlowMove(p : Planet) : TCompFraction;
begin
    if OuterPlanet(p) then
        result := MeanRevolution(p)
    else
        result := MeanParallax();
end;

function TFivePlanet.SlowStep(p : Planet) : TCompFraction;
Var
    d : Integer;
    ArcLen : TCompFraction;
    ArcInZodiac : TCompFraction;

    ZodiacNum : Integer;
    Remainder : Integer;

    Factor : Integer;
    backStep : Boolean;
    Accum : Integer;
begin
    ArcLen := GetMeanSlowMove(p);
    ArcLenFromApogee(p, ArcLen);
    d := ArcLen.data[1] * ArcLen.Wheel[2] + ArcLen.data[2];
    ZodiacNum := d div 135;
    Remainder := d mod 135;

    ArcInZodiac := TCompFraction.Create(ArcLen);
    ArcInZodiac.data[1] := 0;
    ArcInZodiac.data[2] := Remainder;

    Factor := GetSlowStepFactor(p, ZodiacNum);
    if Factor < 0 then
    begin
        backStep := True;
        Factor := -Factor;
    end
    else
        backStep := False;
    ArcInZodiac.Multiply( Factor );
    ArcInZodiac.Divide( 135 );
    ArcInZodiac.Calibrate();

    Accum := GetSlowStepAccum(p, ZodiacNum);
    result := CreateOuterMoveData(p, [0,Accum,0,0,0,0]);

    if backStep then
        result.Sub(ArcInZodiac)
    else
        result.Add(ArcInZodiac);
end;
//@+node:sonam.20180318221449.1: *3* 遲行定數
//@+at 所謂遲行定數，對於武曜，就是其本身行經過運行快慢脩正之后的準確值，對於
// 文曜，就是地球公轉引起的視差在經過地球運行快慢脩正之后的準確值。
// 
// 其計算方法按照《總滙如意寶》第128頁129頁如下：
// 
//     <遲行定數> := <遲行中數> - <遠地點黃經> +/- <遲步>
// 
// 這裏是加還是減取決於是后半週還是前半週。
// 
//@@c
function TFivePlanet.SlowMove(p : Planet) : TCompFraction;
Var
    ArcLen : TCompFraction;
begin
    result := GetMeanSlowMove(p);

    ArcLen := GetMeanSlowMove(p);
    if (ArcLenFromApogee(p, ArcLen)) then
        result.Add( SlowStep(p) )
    else
        result.Sub( SlowStep(p) );

    result.Calibrate()
end;
//@+node:sonam.20180318221530.1: *3* 疾行步度表
//@+at 疾行步度表刻畫了五曜的疾行在一個週期內的角速度變化
// 
//@@c
// 將后半週的逆序下標映射到前半週順序的下標
function TranxSecondHalf(n : Integer) : Integer;
Const
    TranxTab : array [0..13] of Integer =
        ( 1,0,13,12,11,10,9,8,7,6,5,4,3,2 );
begin
    result := TranxTab[n];
end;

// 這裏將《總滙如意寶》表中檢步序數0對應的表項移動到了表的開頭，對應於0下標，
//  同時將后步的表項加上負號，這樣可以稍微方便判斷前步和后步
function GetFastStepFactor(p : Planet; SiderealNum : Integer;
        SecondHalf : Boolean) : Integer;
Const
    FastStepFactorTab : array [0..4, 0..13] of Integer =
        (   (-53,24,23,23,23,21,21,18,15,11,3,-11,-38,-80),
            (-7,10,10,9,8,6,6,2,1,-3,-6,-9,-11,-16),
            (-3,6,5,5,4,4,2,2,0,-2,-4,-5,-6,-8),
            (-34,16,16,15,14,13,11,7,5,0,-4,-11,-20,-28),
            (-73,25,25,25,24,24,22,22,18,15,8,-6,-30,-99)   );
Var
    i : Integer;
begin
    if SecondHalf then
    begin
        i := TranxSecondHalf(SiderealNum);
        i := i - 1;
        if i < 0 then
            i := 13;
    end
    else
        if SiderealNum = 13 then
            i := 0
        else
            i := SiderealNum + 1; 
    result := FastStepFactorTab[Ord(p)-Ord(Mars)][i];
    if SecondHalf then
        result := -result;
end;

// 這裏將《總滙如意寶》表中檢步序數0對應的表項移動到了表的開頭，對應於0下標
function GetFastStepAccum(p : Planet; SiderealNum : Integer;
        SecondHalf : Boolean) : Integer;
Const
    FastStepAccumTab : array [0..4, 0..13] of Integer =
        (   (0,24,47,70,93,114,135,153,168,179,182,171,133,53),
            (0,10,20,29,37,43,49,51,52,49,43,34,23,7),
            (0,6,11,16,20,24,26,28,28,26,22,17,11,3),
            (0,16,32,47,61,74,85,92,97,97,93,82,62,34),
            (0,25,50,75,99,123,145,167,185,200,208,202,172,73) );
Var
    i : Integer;
begin
    if SecondHalf then
    begin
        i := TranxSecondHalf(SiderealNum);
        i := i - 1;
        if i <0 then
            i := 13;
    end
    else
        i := SiderealNum;
    result := FastStepAccumTab[Ord(p)-Ord(Mars)][i];
end;
//@+node:sonam.20180318221543.1: *3* 疾步
//@+at
// 疾步的計算方法按照《總滙如意寶》的第130頁至133頁如下：
// 
//     <疾步> := <疾行積步> + <曜淨行>
// 
// 其中的疾行積步和曜淨行都曜先通過查疾行步度表得到積步和損益率，其中的
// 
//     <曜凈行> := <損益率> * TGenPlanetData(0,餘數,檢步數的剩餘位)
// 
// 查表的檢步序數的計算方法是：
// 
//     <檢步序數> := <檢步> - <遲行定數>
// 
// 得到的差如果滿半週則減去之，結果的最高位即是檢步序數。
// 
// 在如下情況下需要對損益率和所乘的(0,餘數,檢步數的剩餘位)進行特殊處理：
// 
// * 在<檢步序數>減去半週后，如果大於或等於(0,30,0,0,0,0,0)，則需要
// 給<檢查序數>加上(0,30,0,0,0,0,0)。
// 
// * 在查表得到的是前半週的最后一個損益率和后半週的第一個損益率時，如果檢步序數
// 大於或等於(0,30,0,0,0,0,0)，損益率需要乘以2。這表示在近地點時，曜
// 疾行的速度變化加倍。
// 
//@@c
function TFivePlanet.FastStep(p : Planet; var SecondHalf, backStep : Boolean) : TGenPlanetData;
Var
    d : TGenPlanetData;
    HalfCycle : TGenPlanetData; // 半周
    HalfSidereal : TGenPlanetData;
    SiderealNum : Integer;
    Factor : Integer;
    Accum : Integer;
begin
    if OuterPlanet(p) then
    begin
        d := GenDominator(p, MeanParallax());
        d.Sub(ExtPlanetData(p, SlowMove(p)));
    end
    else
    begin
        d := ExtPlanetData(p, MeanRevolution(p));
        d.Sub(GenDominator(p, SlowMove(p)));
    end;

    HalfCycle := TGenPlanetData.Create(p, [13,30,0,0,0,0,0]);
    HalfSidereal := TGenPlanetData.Create(p, [0,30,0,0,0,0,0]);
    if d.GreatEqual(HalfCycle) then
    begin
        SecondHalf := true;
        d.Sub(HalfCycle);
        if d.GreatEqual(HalfSidereal) then
        begin
            d.Add(HalfSidereal);
            d.Calibrate();
        end
    end
    else
        SecondHalf := false;

    SiderealNum := d.data[1];   // 查表檢步序數
    d.data[1] := 0;

    Factor := GetFastStepFactor(p, SiderealNum, SecondHalf);
    if Factor < 0 then
    begin
        backStep := True;
        Factor := -Factor;
    end
    else
        backStep := False;

    if (SecondHalf and (SiderealNum = 0) or not SecondHalf and (SiderealNum = 13))
        and not d.GreatEqual(HalfSidereal) then
        Factor := Factor * 2;

    d.Multiply(Factor);
    d.Divide(60);
    d.Calibrate();

    Accum := GetFastStepAccum(p, SiderealNum, SecondHalf);
    result := TGenPlanetData.Create(p, [0,Accum,0,0,0,0,0]);
    if backStep then
        result.Sub(d)
    else
        result.Add(d);
    result.Calibrate();
end;
//@+node:sonam.20180318221601.1: *3* 疾行定数
//@+at
// 疾行定數即是該時刻各曜視位置的真黃經，按照《總滙如意寶》第133頁，疾行定數
// 的計算方法是：
// 
//     <疾行定數> := <遲行定數> +/- <疾行盈縮步>
// 
// 這裏的加還是減要看前面計算疾步的過程中，曜是運行到了前半週還是后半週，前半
// 周加，后半週減。
// 
//@@c
function TFivePlanet.FastMove(p: Planet; var Direction : Char) : TGenPlanetData;
Var
    Slow : TCompFraction;
    FStep : TGenPlanetData;
    SecondHalf : Boolean;
    backStep : Boolean;
begin
    Slow := SlowMove(p);
    if OuterPlanet(p) then
        result := ExtPlanetData(p, Slow)
    else
        result := GenDominator(p, Slow);
    FStep := FastStep(p, SecondHalf, backStep);
    if not SecondHalf then
        if not backStep then
            Direction := 'E'
        else
            Direction := 'S'
    else
        if not backStep then
            Direction := 'W'
        else
            Direction := 'N';
    if SecondHalf then
        result.Sub(FStep)
    else
    begin
        result.Add(FStep);
        result.Calibrate();
    end;
end;
//@-others

end.
//@-leo
